{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"https://pymt.readthedocs.io\"><img style=\"float: left\" src=\"../../media/pymt-logo-header-text.png\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run a single model in *pymt*\n",
    "\n",
    "In this example we will run a single model in *pymt*. We will use the water balance\n",
    "model, *Hydrotrend* (written in C). We'll setup a model simulation, run it, and then\n",
    "analyze the output. You will also be given a chance to do your own simulations\n",
    "and explore the HydroTrend model. \n",
    "\n",
    "* [Explore the base-case river simulation](#Exercise-1)\n",
    "* [How does a river system respond to climate change?](#Exercise-2)\n",
    "\n",
    "## HydroTrend\n",
    "\n",
    "HydroTrend is a 2D hydrological water balance and transport model that simulates water\n",
    "discharge and sediment load at a river outlet. You can read more about the model, find\n",
    "references or download the source code at: https://csdms.colorado.edu/wiki/Model:HydroTrend.\n",
    "\n",
    "### River Sediment Supply Modeling\n",
    "\n",
    "This notebook is meant to give you a better understanding of what the model is capable of. In this\n",
    "example we are using a theoretical river basin of ~1990 km<sup>2</sup>, with 1200m of relief and\n",
    "a river length of ~100 km. All parameters that are shown by default once the HydroTrend Model\n",
    "is loaded are based on a present-day, temperate climate. Whereas these runs are not meant to\n",
    "be specific, we are using parameters that are realistic for the [Waiapaoa River][map_of_waiapaoa]\n",
    "in New Zealand. The Waiapaoa River is located on North Island and receives high rain and has\n",
    "erodible soils, so the river sediment loads are exceptionally high. It has been called the\n",
    "*\"dirtiest small river in the world\"*.\n",
    "\n",
    "To learn more about HydroTrend and its approach to sediment supply modeling, you can download\n",
    "this [presentation][hydrotrend_presentation].\n",
    "\n",
    "A more detailed description of applying HydroTrend to the Waipaoa basin, New Zealand has been published in WRR: [hydrotrend_waipaoa_paper]. \n",
    "\n",
    "[map_of_waiapaoa]: https://www.google.com/maps/place/Waipaoa+River/@-38.5099042,177.7668002,71814m/data=!3m1!1e3!4m5!3m4!1s0x6d65def908624859:0x2a00ef6165e1dfa0!8m2!3d-38.5392405!4d177.8843782\n",
    "[hydrotrend_presentation]: https://csdms.colorado.edu/wiki/File:SedimentSupplyModeling02_2013.ppt\n",
    "[hydrotrend_waipaoa_paper]: http://dx.doi.org/10.1029/2006WR005570"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "To start, import numpy and matplotlib.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from tqdm import tqdm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And load the HydroTrend model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymt.models\n",
    "hydrotrend = pymt.models.Hydrotrend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "HydroTrend will now be activated in PyMT.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1\n",
    "\n",
    "### Explore the base-case river simulation\n",
    "\n",
    "For this case study, we will run a simulation for 100 years at daily time-step.\n",
    "This means you run Hydrotrend for 36,500 days total. Later on we will change\n",
    "other input parameters but, for now, we'll just stick with the defaults."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymt.models\n",
    "hydrotrend = pymt.models.Hydrotrend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "config_file, config_folder = hydrotrend.setup(run_duration=100)\n",
    "hydrotrend.initialize(config_file, config_folder)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now see that, indeed, we have set up model simulation that will run for\n",
    "100 years. Notice that HydroTrend has a time step of 1 day so we will be\n",
    "generating daily hydrographs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Start time: {0} {1}\".format(hydrotrend.start_time, hydrotrend.time_units))\n",
    "print(\"Current time: {0} {1}\".format(hydrotrend.time, hydrotrend.time_units))\n",
    "print(\"End time: {0} {1}\".format(hydrotrend.end_time, hydrotrend.time_units))\n",
    "print(\"Time step: {0} {1}\".format(hydrotrend.time_step, hydrotrend.time_units))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we look at HydroTrend's input and output variables, we notice that HydroTrend is\n",
    "slightly unusual in that there are now input variables. There are lots of output\n",
    "variables, though."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hydrotrend.input_var_names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hydrotrend.output_var_names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of this list of output variables, we're going to track water and sediment discharge, sediment\n",
    "concentration and bedload flux. The mapping of these variables to the Standard Names\n",
    "reported by *output_var_names* is given in the following table.\n",
    "\n",
    "| Conventional Name      | Standard Name                                             |\n",
    "| :--------------------- | :-------------------------------------------------------- |\n",
    "| Water discharge        | channel_exit_water__volume_flow_rate                      |\n",
    "| Sediment discharge     | channel_exit_water_sediment~suspended__mass_flow_rate     |\n",
    "| Sediment concentration | channel_exit_water_sediment~suspended__mass_concentration |\n",
    "| Bedload flux           | channel_exit_water_sediment~bedload__mass_flow_rate       |\n",
    "\n",
    "What do you think? Are you able to figure out what the other output variables are?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the ***update*** method to advance *hydrotrend* one time step at a time and the ***get_value***\n",
    "method to retreive values from the model. Feel free to add other variables if you like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_steps = int(hydrotrend.end_time / hydrotrend.time_step)\n",
    "\n",
    "q = np.empty(n_steps)\n",
    "qs = np.empty(n_steps)\n",
    "cs = np.empty(n_steps)\n",
    "qb = np.empty(n_steps)\n",
    "\n",
    "for i in tqdm(range(n_steps)):\n",
    "    hydrotrend.update()\n",
    "    \n",
    "    q[i] = hydrotrend.get_value(\"channel_exit_water__volume_flow_rate\")\n",
    "    qs[i] = hydrotrend.get_value(\"channel_exit_water_sediment~suspended__mass_flow_rate\")\n",
    "    cs[i] = hydrotrend.get_value(\"channel_exit_water_sediment~suspended__mass_concentration\")\n",
    "    qb[i] = hydrotrend.get_value(\"channel_exit_water_sediment~bedload__mass_flow_rate\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now plot the complete hydrograph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time = np.arange(len(qs))\n",
    "plt.plot(time, qs)\n",
    "plt.xlabel(\"Time ({0})\".format(hydrotrend.time_units))\n",
    "plt.ylabel(\"Sediment discharge\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate mean water discharge Q, mean suspended load Qs, mean sediment concentration Cs, and mean bedload Qb.\n",
    "\n",
    "*Note all values are reported as daily averages. What are the units?*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(\n",
    "    (q.mean(), hydrotrend.var_units(\"channel_exit_water__volume_flow_rate\")),\n",
    "    (cs.mean(), hydrotrend.var_units(\"channel_exit_water_sediment~suspended__mass_flow_rate\")),\n",
    "    (qs.mean(), hydrotrend.var_units(\"channel_exit_water_sediment~suspended__mass_concentration\")),\n",
    "    (qb.mean(), hydrotrend.var_units(\"channel_exit_water_sediment~bedload__mass_flow_rate\"))\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hydrotrend.var_units(\"channel_exit_water__volume_flow_rate\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Identify the highest flood event for this simulation. Is this the 50-year flood? Plot the year of Q-data which includes the flood.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flood_day = q.argmax()\n",
    "flood_year = flood_day // 365\n",
    "plt.plot(q[flood_year * 365: (flood_year + 1) * 365])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q.max()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate the mean annual sediment load for this river system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qs_by_year = qs.reshape((-1, 365))\n",
    "qs_annual = qs_by_year.sum(axis=1)\n",
    "plt.plot(qs_annual)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qs_annual.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### How does the sediment yield of this river system compare to the present-day Mississippi River?\n",
    "\n",
    "*To compare the mean annual load to other river systems you will need to calculate its sediment yield. \n",
    "Sediment Yield is defined as sediment load normalized for the river drainage area; \n",
    "so it can be reported in T/km2/yr.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2\n",
    "\n",
    "### How does a river system respond to climate change? A few simple scenarios for the coming century.\n",
    "\n",
    "#### What happens to river discharge, suspended load and bedload if the mean annual temperature in this specific river basin increases by 4 °C over the next 50 years?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymt.models import Hydrotrend\n",
    "hydrotrend = Hydrotrend()\n",
    "\n",
    "config_file, config_folder = hydrotrend.setup(run_duration=100)\n",
    "hydrotrend.initialize(config_file, config_folder)\n",
    "\n",
    "n_steps = int(hydrotrend.end_time / hydrotrend.time_step)\n",
    "\n",
    "q = np.empty(n_steps)\n",
    "qs = np.empty(n_steps)\n",
    "cs = np.empty(n_steps)\n",
    "qb = np.empty(n_steps)\n",
    "\n",
    "for i in tqdm(range(n_steps)):\n",
    "    hydrotrend.update()\n",
    "    \n",
    "    q[i] = hydrotrend.get_value(\"channel_exit_water__volume_flow_rate\")\n",
    "    qs[i] = hydrotrend.get_value(\"channel_exit_water_sediment~suspended__mass_flow_rate\")\n",
    "    cs[i] = hydrotrend.get_value(\"channel_exit_water_sediment~suspended__mass_concentration\")\n",
    "    qb[i] = hydrotrend.get_value(\"channel_exit_water_sediment~bedload__mass_flow_rate\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### How much increase of discharge do you see after 50 years? How is the average suspended load affected? How does the bedload change? What happens to the peak event; look at the maximum discharge event of the last 10 years of the simulation?"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
